---
title: "Analysis of Models with Permutations"
author: "Gina Reynolds and Matthew Winters"
output: pdf_document
---

This file is used in replicating findings for "Foreign Aid Funnel? Placebo-Based Assessment of Aid Flows to  Non-Permanent U.N. Security Council Members" by Evangeline Reynolds (University of Illinois) and Matthew Winters (University of Illinois). 


```{r}
library(xtable)
library(stringr)
getwd()
library(foreign)
library(MASS)
Results=read.dta("09_PermuteModelResults.dta")
names(Results)
dim(Results)
NumPermutation #Fed in from the control file

13*2*(NumPermutation+1)

unique(Results$outcome)
unique(Results$subset)

Results$outcome=str_replace(Results$outcome, "odadisUSA_ln"   , "USA (log)")
Results$outcome=str_replace(Results$outcome, "odadisJPN_ln"   , "Japan (log)")
Results$outcome=str_replace(Results$outcome, "odadisDEU_ln"   , "Germany (log)")
Results$outcome=str_replace(Results$outcome, "odadisGBR_ln"   , "United Kingdom (log)")
Results$outcome=str_replace(Results$outcome, "odadisFRA_ln"   , "France (log)")
Results$outcome=str_replace(Results$outcome, "aiddisIBRD_ln"  , "IBRD Disbursement (log)")
Results$outcome=str_replace(Results$outcome, "aiddisIDA_ln"   , "IDA Disbursement (log)")
Results$outcome=str_replace(Results$outcome, "odadisUN_ln"   ,  "UN (log)")
Results$outcome=str_replace(Results$outcome, "odadisWFP_ln"   , "WFP (log)")
Results$outcome=str_replace(Results$outcome, "odadisUNDP_ln"  , "UNDP (log)")
Results$outcome=str_replace(Results$outcome, "odadisUNICEF_ln", "UNICEF (log)")
Results$outcome=str_replace(Results$outcome, "odadisUNHCR_ln" , "UNHCR (log)")
Results$outcome=str_replace(Results$outcome, "odadisUNTA_ln"  , "UNTA (log)")

Results$Outcome_Subset=paste(Results$outcome, Results$subset)

UniqueOutcomeSubset=unique(Results$Outcome_Subset)
UniqueOutcomeSubset

NumOutcomeSubset=length(UniqueOutcomeSubset)
NumOutcomeSubset
```


```{r}

ObsrvdB=rep(NA,NumOutcomeSubset)
ObsrvdB0=rep(NA,NumOutcomeSubset)
ObsrvdB1=rep(NA,NumOutcomeSubset)
ObsrvdB2=rep(NA,NumOutcomeSubset)

ObsrvdT=rep(NA,NumOutcomeSubset)
ObsrvdT0=rep(NA,NumOutcomeSubset)
ObsrvdT1=rep(NA,NumOutcomeSubset)
ObsrvdT2=rep(NA,NumOutcomeSubset)

ObsrvdP =rep(NA,NumOutcomeSubset)
ObsrvdP0=rep(NA,NumOutcomeSubset)
ObsrvdP1=rep(NA,NumOutcomeSubset)
ObsrvdP2=rep(NA,NumOutcomeSubset)

AdjstP =rep(NA,NumOutcomeSubset)
AdjstP0=rep(NA,NumOutcomeSubset)
AdjstP1=rep(NA,NumOutcomeSubset)
AdjstP2=rep(NA,NumOutcomeSubset)

RejectRate95 =rep(NA,NumOutcomeSubset)
RejectRate95_0=rep(NA,NumOutcomeSubset)
RejectRate95_1=rep(NA,NumOutcomeSubset)
RejectRate95_2=rep(NA,NumOutcomeSubset)

Outcome=c()
RegionSubset=c()
  
for (i in 1:NumOutcomeSubset){
  
  SUBSET=subset(Results, Results$Outcome_Subset==UniqueOutcomeSubset[i])
  Outcome=c(Outcome, as.character(SUBSET$outcome[1]))
  RegionSubset=c(RegionSubset, as.character(SUBSET$subset[1]))
  
  RejectRate95[i] =sum(SUBSET$pv[SUBSET$permutation==1]<.05)/NumPermutation
  RejectRate95_0[i]=sum(SUBSET$pv0[SUBSET$permutation==1]<.05)/NumPermutation
  RejectRate95_1[i]=sum(SUBSET$pv1[SUBSET$permutation==1]<.05)/NumPermutation
  RejectRate95_2[i]=sum(SUBSET$pv2[SUBSET$permutation==1]<.05)/NumPermutation
  
  #print(nrow(SUBSET))
  
  ObsrvdP[i] =SUBSET$pv[SUBSET$permutation==0]
  ObsrvdP0[i]=SUBSET$pv0[SUBSET$permutation==0]
  ObsrvdP1[i]=SUBSET$pv1[SUBSET$permutation==0]
  ObsrvdP2[i]=SUBSET$pv2[SUBSET$permutation==0]
  
  temp=sum(SUBSET$tstat >SUBSET$tstat[SUBSET$permutation==0])/nrow(SUBSET)
  AdjstP[i] = if(temp<.5){2*temp}else{2*(1-temp)}
  temp =sum(SUBSET$tstat0>SUBSET$tstat0[SUBSET$permutation==0])/nrow(SUBSET)
  AdjstP0[i]= if(temp<.5){2*temp}else{2*(1-temp)}
  temp = sum(SUBSET$tstat1>SUBSET$tstat1[SUBSET$permutation==0])/nrow(SUBSET)
  AdjstP1[i]= if(temp<.5){2*temp}else{2*(1-temp)}
  temp = sum(SUBSET$tstat2>SUBSET$tstat2[SUBSET$permutation==0])/nrow(SUBSET)
  AdjstP2[i]= if(temp<.5){2*temp}else{2*(1-temp)}
  
  ObsrvdT[i] = SUBSET$tstat[SUBSET$permutation==0]
  ObsrvdT0[i] = SUBSET$tstat0[SUBSET$permutation==0]
  ObsrvdT1[i] = SUBSET$tstat1[SUBSET$permutation==0]
  ObsrvdT2[i] = SUBSET$tstat2[SUBSET$permutation==0]
  
  ObsrvdB[i] = SUBSET$estimate[SUBSET$permutation==0]
  ObsrvdB0[i] = SUBSET$estimate0[SUBSET$permutation==0]
  ObsrvdB1[i] = SUBSET$estimate1[SUBSET$permutation==0]
  ObsrvdB2[i] = SUBSET$estimate2[SUBSET$permutation==0]
  
  B=c("estimate", "estimate0", "estimate1", "estimate2")
  T=c("tstat", "tstat0", "tstat1", "tstat2")
  P=c("pv", "pv0","pv1","pv2")
  Title=c("All", "Unimportant", "Somewhat Important", "Important")
  
  for (j in c(1,4)){  #for now only print All SC and Important
  
   truehist(SUBSET[,T[j]], 
            xlim=c(-5,9),
            ylim=c(0,.5),
            h=.5,
            axes=FALSE,
        xlab="Distribution of t-statistics produced by UNSC Permutations",
        main=paste(Title[j],"Years - ",Outcome[i], "\n",RegionSubset[i]),
        col="white"
        #,prob=FALSE
        )
   axis(2, las=2)
   axis(1, at=c(-4:4))
   abline(v=SUBSET[,T[j]][SUBSET$permutation==0], lty=4)
   abline(v=median(SUBSET[,T[j]]), lty=5)
   legend("topright", 
          #bg="white",
          lty=c(5,4,2,1), 
          legend=c("Median t-stat of \n permuted UNSCs \n  ", 
                   paste0("t-stat of \n realized UNSCs (", round(SUBSET[,T[j]][SUBSET$permutation==0],2) ,")"), 
                   "theoretical t-stat \n distribution",
                   "\n Kernel density \n resulting from estimation \n with permutation"), bty="n",
                    )
      
 #Theoretical Distribution
      x=seq(-5,5,.01)
      y=dt(seq(-5,5,.01), df = 100)  
      lines(x,y, lty=2)    
   
       #Labeling
      #text(2.5,.4 , "Theoretical t-distribution")
      #arrows(2, .387, .8, .34,length = 0.1, angle = 30)
 
   #Observed Distribution
      Kernel<-density(SUBSET[,T[j]], bw=.375)  #smoothing
      lines(Kernel, lwd=1) 
   
  }
  
}

```

\section{Results compared to Vreeland and Dreher}

```{r}
options(width=140)
options(scipen=999)
DF=round(
          data.frame(
             ObsrvdB , ObsrvdT , ObsrvdP ,        AdjstP ,
             ObsrvdB0, ObsrvdT0, ObsrvdP0,        AdjstP0,
             ObsrvdB1, ObsrvdT1, ObsrvdP1,        AdjstP1,
             ObsrvdB2, ObsrvdT2, ObsrvdP2,        AdjstP2,
             RejectRate95, RejectRate95_0, RejectRate95_1, RejectRate95_2
            )
         ,3)
DF=data.frame(Outcome, RegionSubset, DF)

```



```{r, warning=FALSE}


for (i in c("Global", "Africa")){
for (j in c("All Security Council Years", "Important Security Council Years")){

  TEMP=DF[DF$RegionSubset==i,]
  TempRejections=DF[DF$RegionSubset==i,]
  
if(j=="All Security Council Years"){
TEMP=TEMP[c("ObsrvdB", "ObsrvdT" , "ObsrvdP" , "AdjstP") ]
TempRejections=TempRejections[c( "RejectRate95")] 
}
if(j=="Important Security Council Years"){
TEMP=TEMP[c("ObsrvdB2", "ObsrvdT2" , "ObsrvdP2" , "AdjstP2") ]
TempRejections=TempRejections[c("RejectRate95_2")] 

}

names(TEMP)=c("V & D Estimate", "V & D Realized t-stat", "V & D p-value", "Randomization p-value")
names(TempRejections)=c("Rejection Rate at the 95% Level")
rownames(TEMP)=unique(Outcome)
rownames(TempRejections)=unique(Outcome)

print(TEMP)

print(TempRejections)

  }
}

```


```{r,results='markup'}
DF
```

